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ABSTRACT 

In  conventional  inversion  methods,  imaging  structure  inside  the  earth  requires 
reasonable  background  velocities.  In  this  paper,  velocity  analysis  and  structural 
imaging  are  done  at  the  same  time.  The  medium  is  assumed  to  consist  of  constant- 
velocity  layers  separated  by  arbitrary,  smooth  interfaces.  The  objective  of  the 
inversion  is  to  determine  layer  velocities  and  locations  of  the  interfaces.  The 
velocity  analysis  is  based  on  the  principles  that  the  images  will  be  distorted  when 
erroneous  velocities  are  used.  In  particular,  the  difference  between  the  depths 
computed  by  inversions  from  different  experiments  can  be  a  measure  of  the  error  in 
velocity.  The  formulas  of  sensitivity  to  velocity  error  are  derived  for  some  special 
cases.  Some  computer  implementations  for  both  synthetic  data  and  experimental 
data  are  done. 


1.  INTRODUCTION 

Seismic  data  inversion  process  may  be  described  as  follows  :  given  a  scattered 
field  and  a  background  velocity,  reflector  maps  of  the  earth  can  be  constructed  as 
singular  functions  of  reflectors  (  Bleistein,  1987).  The  accuracy  and  efficiency  of  this 
method  for  a  structural  image  of  the  earth  depends  largely  on  the  complexity  of  the 
medium  and  on  the  quality  and  complexity  of  the  description  of  background  velocity. 
Now,  we  suppose  that  the  medium  is  made  up  of  constant- velocity  layers  separated  by 
arbitrary  smooth  interfaces.  Unfortunately,  we  have  limited  information  from  which 
to  guess  these  velocities. 

The  inversion  process  is,  in  fact,  a  form  of  prestack  depth  migration.  Seismic 
records,  as  input,  contain  information  on  travel  time  of  reflected  waves.  Using  a 
given  velocity  model,  we  map  an  image  from  the  time  domain  into  the  depth  domain. 
Use  of  incorrect  background  velocities  generally  results  in  a  distorted  image  of  the 
structure  .  The  degree  of  distortion  depends  on  the  magnitude  of  the  velocity  error 
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and  on  the  length  scale  over  which  signals  are  propagated  with  the  initial  velocity. 
When  the  velocity  is  correct,  the  image  location  should  be  at  the  same  depth,  re¬ 
gardless  of  the  positions  of  sources  and  receivers,  for  the  data  used  in  the  inversion. 
Otherwise,  for  the  incorrect  background  velocity,  the  structural  images  for  different 
source-receiver  offsets  will  differ  from  one  another.  Such  a  deviation  can  help  us  cor¬ 
rect  the  velocity.  After  the  velocity  is  corrected,  the  interface( reflector)  is  determined 
by  picking  amplitudes  from  the  output  image.  To  guarantee  a  successful  ray  tracing 
which  is  used  in  the  inversion  code,  the  modeled  interface  must  be  made  smooth.  This 
process  is  then  repeated  for  determining  the  other  velocities  and  interface  shapes  for 
successively  deeper  layers. 

The  procedure  is  as  follows: 

( 1 )  given  an  initial  guess  for  velocity  model 

(2)  velocity  analysis  on  a  fixed  output  trace  to  correct  velocity 

(3)  using  the  correct  velocity  to  image  an  interface  on  the  output 

(4)  smoothing  the  interface 

(5)  repeating  steps  (1)  to  (4)  for  the  next  layer 

2.  MATHEMATICAL  PRINCIPLES 

We  consider  the  two  dimensional  situation.  We  shall  denote  by  A'  a  2-D  vector, 
X  =  (x,  z)  .  Let  A',  =  Aj(£)  be  source  positions  and  Xr  —  A’r(£)  be  receiver  positions 
located  on  the  datum  surface  L ,  where  £  is  a  position  parameter  on  the  L.  For  any 
point  below  the  surface,  r(Aa,  A)  or  r(A,  Ar),  respectively,  denote  traveltimes  from 
A',  to  A,  or  A’  to  A'P. 

Suppose  we  know  the  total  reflection  travetime  T(£).  Any  reflection  point  X  = 
(x,  z)  must  satisfy 

T(Xa,X)  +  r(A,  Ar)  =  T(£).  (1) 

For  each  the  solution  of  equation  (1)  is  a  curve.  When  f  varies,  we  obtain  a  family 
of  curves. 

Theorem  1.  For  any  velocity  function,  the  envelope  of  a  solution  family  of  equa¬ 
tion  (1)  is  just  the  reflector,  z  =  /(x),  resulting  in  traveltime  T(£). 

Proof.  Differentiate  equation  (1)  with  respect  to  x;  then 

.dr(Aj,  A)  dr(X,  Ar),  3r(A„A)  dr(X,Xr)  dz 

1  dx  dx  1  1  dz  dz  ]dx  ' 

i.e., 

[V*  r(Xs,X)  +  V*  r(A,  Ar)]  ~  =  0. 

This  is  just  a  statement  of  the  law  of  reflection.  Hence  the  curve  z  =  f(x)  is  tangent 
to  the  curve  family  in  equation  (1)  at  every  intersection  point  of  the  family  of  curves 
in  (1)  with  the  reflector.  # 
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By  the  definition  of  an  envelope,  z  =  f(x)  must  satisfy 

t(X„X)  +  t(X,Xt)  =  T(Z).  (3) 

dr(X„X)  dr(X,Xr)  _  dT 

Thus,  for  each  £,  we  can  determine  the  position  of  reflection  point,  (x,  z),  by  equations 
(3)  and  (4).  In  general,  equations  (3)  and  (4)  is  hard  to  solve.  We  only  conseder  some 
special  cases.  Suppose  that  the  medium  velocity  is  a  constant,  c,  and  the  datum 
surface  L  is  the  x-axis.  Then 


A',  =  (xa,  0), 
where 


A'r  =  (xr,0), 


r(X,,X)  =  p,/c , 


r(X,  A'r)  =  pr/c, 


p,  =  ^ /(xa  -  x)2  +  z2,  pr  =  \J(xr  -  x)2  + 

For  this  case,  equation  (3)  and  (4)  are  simplied  to 

p,  +  pr  =  cT(  0, 

dps  ,  dpr  _  r,(  , 
di  +  d£ 

If  we  fix  the  horizontal  coordinate,  x,  of  the  reflection  point,  then  z  and  £  can  be 
considered  as  functions  of  the  velocity  c.  Differentiating  equation  (5)  with  respect  to 
c, 

dp3  dPr  dz  dp,  dpr]dZ  _  ^T>(n,T(n 

1  a7  +  a7>*  + 1  a{  +  a7]*  = cT  (f)  +  T(0 


(5) 


(«3) 


And  using  equation  (6),  we  have 


dz 


(z/p,  +  z/pT)  —  =  T(Z)  =  {p,  +  pr)/c. 


Solving  for  dz/dc ,  then, 


dz  p,pr 


dc  cz 

Introduce  angles  9  and  0  as  in  Figure  1.  Then, 


>  0. 


(7) 


p,  =  z/  cos (9  -  0), 
pr  =  z /  cos(0  +  <t>), 

_ f - .  (8) 

dc  c cos(<?  -  4>)  cos (9  +  4>) 

From  (7),  it  fallows  that  the  imaged  depth  coordinate  of  lefiecuon  for  u.\l-u  x 
(i.e, fixed  trace  location)  is  erroreous,  when  an  incorrect  velocity  is  used.  Moreover, 
the  deviation  is  positive  when  c  is  bigger  than  the  true  velocity  (dc  >  0),  and  negative 
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Fig.  1.  Sketch  of  ray  paths. 


when  c  is  smaller  than  the  true  velocity  ( dc  <  0).  Let  us  discuss  the  relationship 
between  the  deviation  and  the  position  of  the  source  and  receiver.  Let  c*  denote  the 
true  velocity.  At  c  =  c*,  z  is  independent  of  x,  and  xr.  With  no  loss  of  generality,  we 
suppose  that  the  source  is  to  left  of  the  receiver,  that  is,  0  >  0. 


Common-shot  data 


From  Figure  1, 

x3  —  x  =  —  z  tan(0  —  0).  (9) 

If  we  move  from  one  shot  point  to  another,  i.e.  x,  varies,  but  hold  x  fixed,  then 
differentiating  with  respect  to  x,  in  equation  (9)  yields 


1  = 


—z 


de 


cos 2(9  —  <j>)  dx3  ’ 


or 


d9  _  cos2(0  —  4>) 
dx3  z 

By  differentiating  (8)  with  respect  to  x3 ,  using  equation  (10),  and  noticing  that 

dz 

- —  =  0  at  c  =  c  , 
ox. 


(10) 


we  find  that 


d2z 


sin  20 


<0,  at  c  =  c*. 


3x.dc  rctvPM)  4. 

J  \  ‘  T  / 

This  tells  us  that  the  deviation  decreases  as  the  source  moves  right  (Ax,  >  0). 


1 1 : 
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Now,  we  will  obtain  error  estimates  in  velocity  analysis,  which  can  be  derived  from 
(11).  From  figure  1, 


sin  28  _  cos(#  +  0)  _  cos(0  —  0) 
XT  X>  pa  Pr 


Then  (11)  becomes 


d2z 
dx3  dc 


(x3  Xr )  pT 
CZ  Pa 


at  c  =  cm. 


(lb) 


Hence, 

—  ~  .  +  8,2  Z  ■  I  _  .  (c  -  c')  =  iX'-Xr)Pr(c-C’) 

dx3  ~  dx3  c  c  dxs  dc  c  c  c*  z  ps 


(14) 


Suppose  that  we  have  two  shots  x.,, ,  xJ3  and  x,,  <  xS2.  Then  the  difference  in  imaged 
depths  between  the  two  shots  can  be  given  by 


z(xsj  )  -  z(x3l  )  »  dZjX°)  (xS2  -  xsi  )  «  (x30  -  Xr0  )  (xJ2  -  Xat  )  ^  C  ^  P" 


dx3 

where  xSQ  =  (x,2  4-  x„)/2.  Thus, 


c*  2  Pa 


( c  -  £)  _  (2(xJ2)  -  2(xj,,))2  P,  _  A2  2  pa 

c*  ~  (xJ0  -  Xro)(xJJ  -X3l)pr  A Ia(xS0  -  Xro)  Pr' 


(15) 


Obviously,  the  quotient  pa/ pr  is  greater  than  1  for  the  negative  dip  angle,  and  smaller 
than  1  for  the  positive  dip  angle. 

Relationship  (15)  shows  us  the  factors  that  govern  the  accuracy  of  velocity  analy¬ 
sis,  (c  —  c* ) /c* ,  under  the  assumption  that  medium  velocity  is  constant.  The  accuracy 
of  velocity  analysis  is  best  for  large  source-to-receiver  offset,  well-separated  shot  points, 
and  shallower  target.  Interestingly,  it  is  better  also  for  reflectors  with  positive  dip(.i.e., 
receivers  located  in  the  downdip  direction  relative  to  the  shot  point). 


Common-offset  data 

Let  h  be  half  the  offset.  Then, 

2  h  =  z  (tan(0  —  0)  +  tan(0  +  0)). 

Similar  to  the  deduction  of  (10),  we  find  that 

d 2  2  _  2  sin  26  _  , 

dhdc  c(cos2(0  +  0)  -I-  cos2(0  -  0))  > 

Furthermore,  from  (12),  (16)  becomes 

d2z  2  h  2  psPr  . 

ol  3  2i  ^  i 

dhdc  cz  p\  +  pi. 


(16) 


;if) 
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where  h0  =  (h2  +  hi)/2.  Thus, 

(c-c‘)  «  (~(/l2)  ~  z(hi))zfi  +  P*  _  Azz  P2s+P2r  (19) 

c*  ~  2h0(h2-hi)  2pspr  2Ahh0  2p,pr 

The  quotient  {p2  +  p2)/2papr  is  greater  than  1  for  any  given  dip  angle. 

The  relationship  (19)  shows  us  that  the  accuracy  oj  velocity  analysis  deteriorates 
with  increasing  reflector  depth  and  dip,  and  is  best  when  the  two  offsets  and  greatly 
different  from  one  another. 

Note.  From  (19),  we  can  conclude  that  any  error  of  velocity  and  the  difference  of 
the  offsets  result  in  nonzero  deviation  A z.  More  precisely, 

A  /  -wt  l  sd2z(h0)  , 

U-. 

Therefore,  we  define  the  quantity  d2z/dhdc  as  the  sensitivity  to  the  velocity  error. 


Multiple-layer  case 

If  the  medium  is  made  up  of  more  than  one  layer,  the  expression  error  estimate 
must  be  modified.  We  consider  only  a  simple  model,  consisting  of  twro  horizontal 
layers,  and  consider  the  common-offset  situationm  as  in  Figure  2.  Differentiating 
equation  (3)  with  respect  to  C2,  and  using  (4),  we  have 

(dr{X„X)  dr(Xr,X)  dz  _  ,dr(X„X)  dr(Xr,X) 

(  dz  dz  }dc2  1  dc2  dc2  ''  ’ 

For  the  simplicity,  assume  that  Ci  is  equal  to  the  true  value  of  C2.  Then,  6 1  =  d2,  and 


t(Xs,X)  =  r(Xr,X)  =  z/cj  coS#i, 

dr(Xa,X)  _  dr(Xr,  X)  __  cos#i 
dz  dz  ci 

dr(X„X)  ^  dr(Xr,X)  =  -d2 

dc2  dc2  cl  cos 
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From  this  and  (20),  we  find  that 


dz 


dc2  ct  COS#i 


2  ' 


(21) 


Furthermore,  the  sensitivity  to  the  velocity  error  is  given  by 

d2z  .  d2  2tan#i 


dh  dc2 


ICJ=Ci  ' 


d\  +  d2  C1 


(22) 


compared  to  (16)  (let  6  =  0).  Equation  (22)  shows  us  that  for  multiple  layers,  the 
sensitivity  to  the  velocity  error  in  a  thin  layer  is  reduced  by  the  ratio  of  the  layer 
thickness  to  layer  reflector  depth. 


3.  SPECIAL  TECHNIQUES 


Iteration  for  velocity 

If  we  view  the  depth  difference  A z  as  a  function  of  velocity  c,  then  the  correct 
velocity  c‘  is  the  one  for  which 

Az(c)  =  0. 

Based  on  this  observation,  we  use  a  simple  iteration  to  determine  the  velocity. 

(1)  Choose  initial  velocities  q  and  C2  such  that 

dx  =  Az(q)  <0,  d2  —  A2(c2)  >  0. 

(2)  Compute  a  new  velocity  by  weighting  the  initial  velocities  as  follows: 


C3  =  Ci  d2/ (d2  —  di)  +  c2di/(di  —  d2). 

(3)  If  A 2: ( c3 )  =  0(or  smaller  than  a  given  precision),  the  iteration  will  stop,  with 
c3  the  desired  velocity.  If  A2(c3)  >  0,  update  c3  replacing  c2  in  step(l);  otherwise,  if 
A:(c3)  <  0,  update  c3  replacing  q. 

Smoothing  the  interface 

Ray  tracing  in  the  inversion  code  is  stable  when  the  description  of  the  interface 
has  second-order  smoothness.  Consequently,  smoothing  of  the  interface  before  ray 
tracing  is  desirable. 

Let  z  =  f(x)  be  any  continuous  function.  We  solve  for  a  smooth  function  g(x) 
that  approximates  f(x)  through  the  requirement 

J (/(*)  -  g{x))2dx  +  a  j (~)2dx  =  min,  (23) 
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where  a  >  0  is  called  the  smoothing  parameter.  The  larger  the  value  of  a,  the 
smoother  will  be  g( x)  . 

By  calculus  of  variations  we  can  change  (23)  into  a  differential  equation  for  <j{x). 
For  any  positive  number  A  and  any  smooth  function  witn  zero  boundary  condition, 
t j  =  rj(x ),  we  define  a  functional, 

B{\,n)  =  J [/(*)  -  (g(x)  +  '\t?)]2 dx  +  a  /(^  +  \^-~)2dx. 

Then  (23)  is  equivalent  to  jjl  |a=o=  0.  for  any  tj.  That  is, 

J (ff(x)  ~  f{x)}rj(x)  dx  +  a  J  dx  =  0. 

Using  integration  by  parts,  we  have 

J [$(*)  -  /(*)  +  a  — ] tj(x)  dx  =  0. 

Since  tj  is  arbitrary,  this  equality  is  satisfied  if  and  if 

Taking  the  Fourier  transform  gives  the  solution  in  the  wavenumber  domain 

G(k)  =  F(k)/(1  +  ak4). 

This  expression  shows  that  high- wavenumber  components  of  f(x)  are  suppressed  in 
the  approximation  g(x) 

4.  COMPUTER  IMPLEMENTATIONS 

To  testify  the  efficiency  of  our  method,  we  do  a  number  of  numerical  common- 
offset  experiments.  The  inversion  code  is  based  on  the  assumption  that  the  medium 
is  two-and-half  dimensional  (Hsu,  1991). 

Example  1:  Modeling  data 

First,  we  take  synthetic  common-offset  data  from  the  layered  model  shown  in 
Figure  3.  The  input  synthetic  data  were  obtained  by  a  common-offset  modeling 
program.  The  synthetic  data  are  generated  for  five  gathers  with  common  offsets 
100  m,  300  m,  500  m,  700  m,  and  900  m,  with  shots  and  receivers  on  the  horizontal  top 
surface.  The  first  shot  is  at  x  =  100  m,  and  the  shot  point  spacing  is  20  m.  Each  offset 
uses  100  shots  and  receivers.  The  sampling  interval  is  4  ms,  and  the  total  reflection 
time  is  2  sec.  The  inversion  output  spans  the  ranges  200  m  to  1800  m  in  x,  and  0  to 
3000  m  in  z.  Velocity  analysis  is  done  through  the  third  layer  as  in  Figure  6.  After 
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obtaining  the  velocity  in  one  layer,  we  pick  an  interface  from  the  inversion  output 
with  this  approximate  velocity,  then  smooth  it.  The  process  is  repeated  recursively 
through  all  the  layers.  The  results  are  as  follows  : 

(1)  In  the  first  layer,  the  iterative  values  of  Ci  are  1500,  3000,  2013;  the  true  value 
is  2000. 

(2)  In  the  second  layer,  the  iterative  values  of  C2  are  2013,  3500,  2863,  3022;  the 
true  value  is  3000. 

(3)  In  the  third  layer,  the  iterative  values  of  c3  are  3022,  4500,  4007;  the  true  value 
is  4000. 

Us’..g  the  model  consisting  of  these  approximate  velocities  and  interfaces,  we 
obtain  an  inversion  output  which  is  very  closed  to  the  inversion  output  using  the 
exact  model  consisting  of  the  accurate  velocities  and  interfaces.  (See  Figure  5.) 

Now  we  test  the  sensitivity  to  velocity  error.  The  theoretical  values  are  computed 
from  the  equations  (19)  and  (22).  In  all  layers,  take  Ac  =  500  m/s,  A h  =  400  m, 
and  hQ  —  250  m.  In  the  first  layer,  the  theoretical  A:  is  100  m,  the  measured  value 
is  130  m;  in  the  second  layer,  the  theoretical  A z  is  17  m,  the  measured  value  is  30  m; 
in  the  third  layer,  the  theoretical  A^  is  12.5m,  the  measured  value  is  10  m.  The 
least  measurable  Az  (the  sample  spacing)  is  10  m  and  this  limits  the  accuracy  of 
velocity  analysis.  The  errors  between  the  theoretical  and  measured  A:  are  1/4  to  1/2 
wavelength  at  the  dominant  frequency. 

Example  2:  Marathon  data 

The  input  data  is  from  a  physical  experiment.  The  real  medium  can  be  approxi¬ 
mated  to  a  two-and-half  dimensional  model.  The  data  were  296  shots,  each  shot  with 
48  receivers.  The  shot  point  spacing  is  80  ft,  and  the  receiver  interval  is  80  ft.  We 
sorted  the  data  into  five  gathers  of  common-offset  with  offsets  880  ft,  1680  ft,  2480  ft, 
3280  ft,  and  4080  ft.  The  first  shot  point  is  at  x  —  0.  For  each  offset,  there  are  256 
shots  and  receivers.  The  sampling  interval  is  4  ms;  the  total  time  is  2  s.  The  inversion 
output  spans  the  ranges  x  from  200  to  24000  ft  and  z  from  0  to  12000  ft.  Velocity 
analysis  is  done  through  the  fourth  layer.  The  results  are  follows: 

(1)  In  the  first  layer,  the  iterative  values  of  C\  are  8000,  13000,  10857,  11714:  the 
true  value  is  11750. 

(2)  In  the  second  layer,  the  iterative  values  of  c,  are  11714,  17000,  15679;  the  true 
value  is  15750. 

(3)  In  the  third  layer,  the  iterative  values  of  Cj  are  15679,  24000,  21920;  the  true 
value  is  22410. 

(4)  In  the  fourth  layer,  the  iterative  values  of  cj  are  21920,  13000,  15973;  the  true 
value  is  15750. 


9 


Liu  and  Bleistein 


Velocity  Analysis 


5.  CONCLUSION 

In  section  4,  we  did  computer  implementations  for  both  synthetic  data  and  exper¬ 
imental  data.  The  numerical  results  show  that  our  method  can  obtain  high  accu  acy 
in  the  estimate  of  velocity  and  the  imaging  of  interface.  Moreover,  as  a  practical 
inversion  method,  there  exist  some  questions  to  be  stressed  on. 

Model  limitations 

The  present  inversion  code  requires  that  the  medium  be  made  up  of  constant- 
velocity  layers  separated  by  smooth  interfaces.  But,  in  actual  subsurface,  interfaces 
may  touch  each  other  or  terminate  abruptly,  and  velocities  may  vary  laterally.  Usu¬ 
ally,  we  separate  these  interfaces  on  purpose  in  order  to  guaraniee  a  successful  inver¬ 
sion  and  assume  velocities  are  constants.  However,  it  will  produce  the  model  error. 
To  solve  this  problem  completely,  a  new  inversion  code  for  more  general  models  needs 
to  be  devised. 

Selection  of  the  output  trace 

In  theory,  one  output  trace  with  a  fixed  horizontal  coordinate  is  sufficient  to  de¬ 
termine  the  velocity.  However,  if  the  error  of  model  cannot  be  ignored,  such  estimate 
of  velocity  may  be  unstable.  A  better  way  is  to  select  more  one  output  traces.  Then 
take  the  average  value  of  these  velocities  from  the  different  traces. 
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Fig.  7.  Model  tank  data  for  two  offsets. 
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Fig.  8.  Inversion  outputs  for  the  data  in  Figure  7.  The  input  model  was  obtained 
from  velocity  analysis.  Five  offsets  were  used  in  the  stacking. 
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Fig.  9.  Velocity  analysis  on  the  model  tank  data.  The  true  velocity:  c\ 
11750,  c2  =  15750,  C3  —  22410,  c4  =  15750. 
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